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List of Symbols 


The item in parentheses indicates the equation number or section where the 
symbol first appears. 

a,b,c coefficient matrices of the nonconservative Euler equations, (App. C) 

c speed of sound, (9); chord length, (Sect. VI) 

e energy per unit volume, (1) 

f,g,h flux vectors in Cartesian coordinates, (1) 

i,j indices representing £,n,s directions, (Mb) 

i,j,k representative curvilinear coordinates, (9) 

2,,ra indices representing elements of the flux Jacobian matrices, (6b) 

n normal to cell face, (2c); time level, (5) 

p static pressure, (1) 

q conserved dependent variables in Cartesian coordinates, (1) 

nonconserved dependent variables, (C3) 

t time, (1) 

u,v,w velocity components in Cartesian coordinates, (1) 

Ug boundary layer friction velocity, (Sect. IV) 

x,y,z Cartesian coordinates, (1) 

x representative element of the dependent variable vector, ( 23 ) 

A ,B ,C flux Jacobian matrices, (6b) 

A representative element of the coefficient matrices, (21) 

CP pressure coefficient with respect to freestream conditions, (Figures) 

F,G,H flux vectors in the transformed equations, (2) 

I,J,K indices for computational grid lines, (Sect. VI) 

H shape factor, (Sect. IV) 

I identity matrix, (8b) 
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J Jacobian matrix of the transformation from Cartesian to curvilinear 

coordinates, (2) 

K representative flux vector (App. C) 

K representative flux Jacobian matrix (App. C) 

M Mach number, (31); 3Q/3q, (C4) 

P eigenvectors of nonconservative Euler equations, ( C 1 1 ) 

Q conserved dependent variables in the transformed equations, (2) 

R residual, (13) 

Re- Reynolds number based on mean aerodynamic chord, (Sect. VI) 

5 cell surface area, (2c) 

T transpose, (1); eigenvectors of the conservative Euler equations, (CIO) 

U,V,W contravariant velocity components, (2) 

X representative unknown dependent variable vector, (21) 

X I ,X 2 unknown dependent variable vector in the two pass algorithm, (14) 

a angle of attack, (Figures) 

Y ratio of specific heats, (1) 

6 spatial difference operation, (3) 

6 boundary layer displacement thickness, (Sect. IV) 

5,n»C curvilinear coordinates, (2) 

0 contravariant velocity, k x u+k y v+k z w, (C7) 

k representative coefficient matrix for the nonconservative Euler 

equation, (App. C) 

\ eigenvalues, (9) 

v cell volume, (2a) 

p density, (1) 

t time in transformed equations, (2) 

4> (u 2 + v 2 + w 2 ), (C7) 


i 
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A difference operator (usually first order), (2a) 

A diagonal matrix of eigenvalues, (Cl 5) 

V del operator, (2a) 

( )„ freestream conditions, (31) 

C) ( )/ | Vk | , (C17) 

( ) e boundary layer edge conditions, (Sect, IV) 

( ) vector quantity 

( )- terms associated with positive or negative eigenvalues, (10c) 
(~) equivalent incompressible boundary layer terms, (Sect. IV) 


I . Introduction 


During the mid 1970's a computer program for the solution of the 
Euler equations was developed by researchers at Sverdrup Technology, Inc. 
for transonic and supersonic applications. This program, known as ARO- 
1,(15,16) was an explicit, finite volume formulation based on McCormack's 
predictor-corrector algorithm. From this code has evolved a long line of 
Euler equation solvers whose current members now share only the finite 
volume formulation with the prototype. One of these current programs is 
the subject of the following presentation. 

The code of interest here goes by the (rather colorful) name BROWN 
MULE. It is a three-dimensional, time accurate, finite volume Euler 
solver (as its progenitors) with extensive refinements. ^ » 2 »3) The solu- 
tion scheme is an implicit, upwind, split-flux vector formulation in which 
the flux vectors are divided into subvectors based on the signs of the 
eigenvalues. Prudent approximate factorization then leaves only a system 
of block bidiagonal equations to be solved which is readily accomplished. 
The upwind scheme used here requires no artificial dissipation and is 
conditionally stable in three-dimensions. Furthermore, because the pres- 
ent interest is in steady state solutions, local time stepping can be used 
thereby improving the convergence rate. The resulting program is easy to 
use, requiring only a minimum of input variables and parameters, and 
achieves engineering quality answers in reasonably short machine time. A 
disadvantage is the relatively large memory needed to store the flow 
variables and split flux vectors, although with the latest generation of 



supercomputers this is not a serious limitation. Further history and 
background concerning BROWN MULE and its family may be found in References 
1 and 3 . 

As part of the present research, a version of BROWN MULE was devel- 
oped to include calculation of viscous effects, including effects of 
moderate flow separation. This is accomplished through an unique inverse 
integral boundary layer solution which uses an analytical description of 
the velocity profile. This profile and other supporting relationships 
are based on curve fits to experimentally determined compressible turbu- 
lent boundary layers including separated layers. The viscous effects are 
imposed upon the flow through a surface source model which is simply imple- 
mented as a modified solid wall boundary condition. These viscous calcu- 
lations are quite efficient and impose negligible time and storage 
penalties on the overall program so that the capability of the program is 
significantly enhanced at little cost. 

This paper explains the theoretical and numerical bases of the pro- 
gram with emphasis on the logic behind the equation development. In 
addition, the program is fully detailed so that a user can quickly become 
familiar with its operation. 

Because this code is intended for computation of complex flow fields, 
an application to transonic flow past a wing/body configuration represen- 
tative of a modern wide body turbofan transport is also presented. A 
companion paper describes in detail the explicit vectorization of this 
program on the NASA Langley Research Center VPS-32. 
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II. Formulation of the Numerical Procedure 


Governing Equations 

The governing equations for inviscid flow, the Euler equations, take 
the form 


ia + ll + ls + i*} . o 

3t 3x 3y 3z 


( 1 ) 


where 


q - [p, pu, pv, pw, e] T 


f - [pu, pu 2 + p, puv, puw, u(e+p)] T 


[pv, puv, pv 2 + p, pvw, v(e+p)] T 


[pw, puw, puv, pw 2 + p, w (e+p)] T 


„ . i r ^ r 2 . 2. 2m 

P * ». i - 1 u e - -jjp u v ♦ w jj 


This form of the equations (other forms are possible with algebraic ma- 
nipulation) is said to be in strong conservation form ^ 1 ^ because the funda- 
mental properties conserved in nature - mass (p), momentum (pu, pv, pw) and 
energy (e) - always appear as distinct entities in the equations. The 
transformation ^ (outlined in Appendix A) using curvilinear coordinates, 
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5 - 5 (x t y,z,i) 


n * n(x,y,z,T) 


C - c(x,y,z,T) 


T = T 


leads to 


3Q + 3F + 3G + 3H 

9t 35 3n 3; 


( 2 ) 


where 

Q = J [p , pu, pv, pw, e] T 

F = J [pU, puU + 5 x p, pvU + 5 y P, PWU + 5 Z P* U(e + p)] T 

G = J [pv, puV + n x p, pW + n y p, pwV + n z p, V(e + p)] T 

H ** J [pW, puW + c x p, pvW + c y p + pwW + c z p, W(e + p)] T 

J “ x 5^ y n z C ” z n y C^ " y 5^ x n z C ” z n x c^ + z c^ x n y c ~ y n x ^ 


U - 5 x u + 5 y v + 5 z w + 5 t 


V = n x u + riyV + n z w + n t 
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w = s x u + c y v + e z w + s t 


5 x - J "' ( Vc ‘ Vc> 


"x ' J '’ [ \h ' Vc> 


5 y ' J "' ( Vt ' \\ ] 

5 z ' J ‘' ( Vc ' Vs* 


' \ 5 x ' y , 5 y - z t E x 


j " 1 (x_z - z.xj 


rc i c‘ 
C ^ y 5 " W 


n„ - J 1 (x r y. - y x_) 


\ - “Vx ' Vy ” Vz 


x, = J 1 (y z 
x n 

S ' J "’ ( V« 
5 z ’ J '’ ( Vn 


Vn> 


Vs 1 


Vn 5 


Discretization 


? t = “Vx - Vy - Vz 


To discretize (2), let us integrate equation (2) over the volume of a 
computational cell. 


f r 3Q A 3F A 3G ^ 3H-, . 

J ^ + 3? + 3^ + 3? )dV 


(2a) 


The mean value theorem of calculus permits us to write 


f 3Q , r AQ'| 

J 3t dv = 
v 

- AQ^ 


( 2 b) 


where (^~) is an effective average value of within the cell 

volume v. 

The flux terms combine to give /V* Fdv, where 7 is the divergence operator 


(in the curvilinear coordinate frame) and F is a vector with components F, 
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G and H. 

The divergence theorem then permits us to write 

fV’Fd = [n*FdS (2c) 

J S J S 

where S is the cell surface and n is the unit outward normal. 

The flux terms - F, G, H - transport fluid properties across C,n,c faces 
of a cell, respectively, so that (2c) becomes 

j s UdS - (F out - FjiniC * (G out - G ln ) 4C 4 C * (H out - HjiSJn (2d) 

where effective values of F, G and H are used on the cell 
faces. Using (2b) and (2d) in (2a), letting 6 = ( ) Qut - ( ) in , realizing 
Av = A£AnAc, and dividing through by Av gives 


AQ + 6F + 6G + 6H 
At A5 Ar> Ac 


( 3 ) 


Equation (3) might have been obtained from equation (2) by simply approxi- 
mating the partial derivatives in (2) with suitable differences. However, 
the procedure [(2a) through (2d)] here shows more clearly the finite volume 
nature of equation (3) and, therefore, of the entire problem. The volume 
approach is, in general, a more physically meaningful approach and is 
recommended whenever the option for this approach arises. It is now neces- 
sary to precisely define the individual operators and terms in equation 
(3) and to formulate the precise numerical problem to be solved. The 6 
operator is the central difference operator notation and means, e.g. 

6 i F = F i+1/2 ~ F i-1/2 (4a) 
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To interpret this relation we realize that the cell center is at i, so 
that i + 1/2 and i - 1/2 signify cell faces, and thus the flux vectors 
(F, G, H) are evaluated at cell faces as dictated by the finite volume 
approach. The computational coordinates ((■, n, ?) correspond directly to 
the cell indices (i, j, k) so that A? = An = A? = 1 and (3) can be written 


AQ = -At (fijF + 6 j G + 5 k H) (4b) 

We will be advancing the solution in time. The most recently computed 
time step is the n step, or level. Equation (4b) might then be written as 


AQ n «■ -At (6 i F n+1 + <5jG n+1 + <5 k H n+1 ) (5) 

with AQ n - Q n+1 - Q n 

The unknowns are really Q n+1 , but it is convenient to work with AQ n as the 
unknowns because eventually, as we get closer to the final solution, AQ n 
will go to zero. This fact will be used to improve our computational 
efficiency. 

Linearization 

Now, the flux vectors are explicit functions of Q, so that we can 

write^'^*®) 

F n+1 
G n+1 
H n+1 


F n ♦ (f|) n (Q n+1 
G n + ( |Ojn (Q n + 1 

H n ♦ (f ) n (Q n+1 


Q n ) 

Q n ) 

Q n ) 


+ • 




+ • • • 


F n + (|^) n AQ n + - 


3Q' 

c n * (|2) n 4Q n 


+ • • • (6a) 


Hn + AQ " + 
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The flux vectors are nonlinear in the dependent variables, Q or &Q. For a 


manageable numerical solution they must be linearized, which is accom- 
plished by dropping terms of order (AQ n )^ and higher. The derivatives 

such as 3F/3Q, which appear in the expressions for the flux vectors, form 
matrices known as the Jacobian matrices for each flux vector. Because 
equation (1) or (5) actually represents 5 distinct equations, there are 5 
elements to each of Q, F, G and H and the flux Jacobian matrices therefore 
each contain 25 elements, 


Im 


K_ 

3Q n 


, B 


im 


3G ? 

3Q„ 


, C 


n 

£m 


jsun 

A , 1 - 1, 5, m 
3Q„ 


1, 


'm nn ^m 

With this notation, the linearized flux vectors become, 

_n+1 _n , , n ._n . ,n ._n , ,n .-n . n . rt n . .n t -n 

F 1 =F 1 + A 11 AQ 1 + A 12 AQ 2 + V°3 + A 14 AQ 4 + A 15 AQ 5 


(6b) 



F n + A n AQ n ♦ A n AQ n + 
5 51 1 52 w 2 


+ A 


55 



with similar expressions for G n+ ^ and H n+ ^ . 

Using (6b), and the companion expressions for G and H, in (5) gives 


AQ 


AQ* 


♦ J, (*!>,*...♦ *" 5 4Q 5 n )) ♦ 8j 0j * 

* 5 k H ? * < k( c n M ?*-* c ?5 4Q ? ) 

(7) 

-i,( Sl F 5 n * Sj (A^-* »" 5 AQp * A jG^ * Sj (b",AQ^...*b” 5 Aq") 
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Equations (7) are an implicit set of linear equations for the unknowns 

AQ n . They are implicit in two ways. First, each AQ n depends on all five 

AQ n , giving a system implicit with respect to the dependent variables at a 

given location in the grid. Second, the difference operator introduces 

neighboring values of the unknowns (recall, for example, 

. „n+1 „n+1 „n+1 . , . . 

6 F i = Fj -F i-i thus Eivmg a 

2 2 

spatially implicit system as well. Equations (7) can be algebraically 
rearranged into the following matrix equation. 


1 * 4t ( s i A n* 5 j B n* s k c n) 

• 

At(6 .a? ,+<5 .B? +6 c" ) 

1 i 12 j 12 k 12 ; 

• 

AT ^ 6 i A 52 +l5 j B 52 + ' 5 k C 52^ 


• AT[6 i A 15 +6 j B 15 +6 k C 1 5 f Q 1 


AQ 

|aq 

AQ 


(8a) 


.1+At 




-mKs.F^SjC^S^") (SjFJ* «j0^ k H 2 n ). 


•K^V5*v 5 n )l T 


or, in compact notation 

(I+At<5 . A *+At5 .B *+At6 O ) AQ n = -At( 6 .F n +5 ,G n +6 H n ) (8b) 

1 J K X J K 

The dot indicates that the difference operator acts on the products AAQ n , 
BAQ n and CAQ n . 

Equations (8) are the expression of the Euler equations which are to 
be numerically solved. It might be wise, at this point, to review the 
basic steps leading to (8). We start with the Euler equations in conser- 
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vative differential form in Cartesian coordinates, equation (1). These are 
transformed to curvilinear coordinates, equation (2). The transformed 
equations are discretized by integrating over a cell volume to yield 
equation (3). Details of the discretization are introduced and finally 
the flux vectors are linearized with respect to the dependent variables. 
The result of these last operations is equation (8). Thus we see that (8) 
is a discretized, linearized, finite volume formulation of the Euler equa- 
tions. 

Eigenvalues 

The problem now becomes one of solving the system of equations in (8). 
Various schemes are possible; the scheme used here is centered on the 
notion that information is propagated through a flow field in certain 
preferred, or characteristic directions, and with characteristic veloci- 
ties. (The word ’characteristic', as used here, has a double meaning. 
The first meaning may be taken as 'particular' and is motivated by physi- 
cal arguments, while the second meaning stems from the mathematical prop- 
erties of systems of partial differential equations.) Physically, 
information is propagated via bulk fluid motion and via acoustic waves. 
Acoustic waves travel in all directions between molecules at the local 
sound speed; the molecules, meanwhile, are being convected, on the aver- 
age, at the bulk flow velocity. Thus, the net acoustic wave speed is the 
sum of the sound speed and bulk velocity. The fluid velocity and the net 
acoustic wave speeds are the characteristic velocities. These concepts 
are mathematically implemented through consideration of the eigenvalues of 
the flux Jacobian matrices. 


10 



The eigenvalues of the flux Jacobian matrices are the mathematical 
characteristics of the system represented by (8). For the present problem 
they turn out to be identical to the characteristic velocities with which 
information is propagated in the flow. (The theory of characteristics is 
described, at length, in Ref. 7; the operations required to find the 
eigenvalues of the matrices here are outlined in Appendix C). Each flux 
matrix is a 5 x 5 system and has, therefore, 5 eigenvalues 


A*, A 1 , A* , i * 1, 5 with (5,n,c) referring to (A,B,C), respectively, 
s "H C 

The eigenvalues for the flux Jacobian matrices here are 

A. 1 = X 2 = A^ * k u + k v + k w 
k k k x y z 


- »k * °l Vk l 


(9) 


‘k - l k - 'i 7k i 

where k = £,n,c for F,G,H respectively, c is the speed of sound and 


|Vk| = (k* + k 2 * k 2 z ) U2 . 


Interpretation of the eigenvalues as measures of the velocities at which 
fluid properties and information are propagated through the flow field is 
evident in these expressions. 

Flux Vector Splitting 

To take advantage of the physical significance of the eigenvalues, we 
write the flux vectors as a linear combination of, so-called, subvectors 
which have as coefficients the five eigenvalues of each vector.^ »2»3) 
Hence 
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i = j = 1 X 5 fi J = X 5 f i1 + X C f i1 + A 5 f i3 + + X 5 f i5’ 1 1,5 


(10a) 


and similarily for 
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Since, X, "* = 
k 

to give = 


X i s i i 

n ^ j 


X^ f 

Vil 


and H = f X'jh 
j=1 C J 

3 

X , then the first three terms can be combined together 

li S 

+ X^f^ + and similarly for G and H. (10b) 


A very important sequence of steps is now undertaken. Recall that the 
equations (8) are spatially implicit. This implicit aspect of the equa- 
tions can be removed as follows. The flux vectors F,G,H, by (10b) are 
sums of subvectors based on the eigenvalues. These subvectors can be 
grouped together according to whether their eigenvalue coefficients are 
positive or negative, so that 

F = F + + F", G = G + + G~" and H = H + + H~ (10c) 

where, for example, F + is a subvector made up of the Xf terms which have 
X>0 and F~ is composed of the remaining Xf terms (with X<0). 

A simple physical idea motivates the above decomposition of the flux 
vectors into components with positive and negative eigenvalues. Because 
the eigenvalues are the characteristic velocities of information propaga- 
tion, then their sign indicates in which direction information is carried 
through the flow field. For example, suppose 5 increases to the right. 
Then, if X^ is positive this means it is carrying information to the 
right. Conversely, a negative value of X^ indicates information moving 
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to the left. Now F represents the actual transport of information through 
the flow. Since X ^ can be positive and negative, then F can transport to 
a cell from both the left and the right. A similar argument applies to G 
and H so that, in general, information reaching a cell comes from all 
surrounding cells. The surrounding cells, themselves, contain unknown 
qualities so that the problem is spatially implicit; the solution at each 
cell depends on the solution in all the surrounding cells. However, if we 
decompose the flux vectors into components with positive and negative 
eigenvalues then we can isolate contributions to a cell due to the flux 
from any particular direction. The ability to isolate flux contributions 
is used to develop a numerical algorithm which, at each stage in the 
calculation, involves only flux terms transporting known quantities into 
each cell and thereby removes the spatially implicit nature of the equa- 
tions. 

Returning now to equation development, we observe that each of these 
new subvectors has a Jacobian matrix associated with it, that is 



These new subvectors can also be linearized exactly as the entire vector 
was linearized to yield equation (6b), 

(F + ) n+1 = (F + ) n + (Q n+ ^“ Q n ) and so on for the other subvectors. 

We now take this linearization plus equations (10) and (11) to form an 
equation analogous to (8), 
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Approximate Factorization 


Equation (12) is a block tridiagonal system and is still implicit, 
both spatially and with respect to the dependent variables at a given 
point. To eliminate the spatial implicitness, an approximation to equa- 
tion (12) is formed via approximate factorization,^ 

(i + Ax6 i A* + AxA^B* + AxA^C* ] (i + AxA.A~ + AxAjB~ + At^c") AQ n = -AxR n 

(13) 

where R n = residual = (A.F + + A.F* + A .G + + A 4 G~ + A, H + + A,H~) n 

v i l j j k k 

In this form, terms of order (At) 2 appear on the left hand side which are 
not in the original equation (12). Consequently this formulation is, at 
best, second order accurate in time. We solve equation (13) here with the 
two step scheme 

(I + AxA^A* + ATjB* + AxA^C* ) X 1 = -AxR n 
(I + Ax6 i A~ + AxA.B* + Ax<$ k C“) X 2 - X 1 (14) 


Final Equations 

The algorithm in equation (14) is now illustrated. Consider the first 
equation multiplied out, 

X 1 + AxA (aV) + AxA . (B + X^ ) + AtA (cV) = -AxR n 

i j k 

and write out the finite difference operations using a first order spatial 



difference. 


x 'i,j.k* 4t(AV) i.J.k- 4 *< #V >i„,j.k * 4 *< BV >i,j.k ' 4 * (BV) i„r,,, 


‘^‘A.J.k-^’Ji.j.k-, -- 4 <j. 


(15) 


The flux matrices (A + ,B + ,C + ) are functions of the cell metrics and the 
dependent variables. The metrics are evaluated at the cell faces indi- 
cated by the subscripts appearing in equation (15). The dependent vari- 
ables to be used in evaluating the flux matrices are determined based on 
the facts that (1) this is a cell centered, finite volume scheme so that 
the dependent variables are taken as constant within a cell, and (2) the 
matrices (A + ,B + ,C + ) involve only positive eigenvalues and represent trans- 
port in the positive (5,n,?) directions only. For example, considering 
the C flux through cell (i,j,k), the flux in comes from cell (i-1,j,k) 
while the flux out is from cell (i,j,k), which means we use 

Q 1 ? . . to evaluate (A + X^). . . and Q?_. . . to evaluate (A + X^)._. . . . 

1 , J , K 1 I »J ^ I ij 


To show this we will express the flux matrices as functions of the depend- 
ent variables whose indices indicate the locations at which they are 
evaluated. Equation (15) now becomes 

X i,j,k + Ax ^ A (Q i,j,k^i,j,k X i,j,k " Ax ^ A (Q i-1 ,j,k^i-1,j,k X i-1,j,k 
+ Ax ^ B (Q i,j,k^i,j,k X i,j,k - Ax ^ B (Q iJ-1,k )] i,j-1,k X i,j-1,k 
+ Ax ^ C (Q i,j,k^i,j,k X i,j,k - Ax ^ C (Q i,j,k-1^i,j,k-1 X i,j,k-1 


-AtR 


i.J.k 


( 16 ) 
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or rearranging, 


U ♦ Ax[A + (Q^ j>k ) ♦ B + (Q^ Jfk ) ♦ C + (Q^ J>k )] if Jik }xl fJfk - 

+ At ^ A (Q i-1 f j,k^i“t,j,k X i-1,j,k + At ^ B (Q iJ-l,k^i f J-l.k X i f j*l,k 


+ At[c (Qj_ , j,k-1 X i,j,k-1 


(17) 


Examining equation (17) reveals that if we start at the lowest index 

boundaries and advance in the direction of increasing (i,j,k) consistent 

with the positive eigenvalues, then the right side of (17) is always known 

and therefore the unknowns, xl . , can be solved for directly. (The 

i * J » k 

boundary conditions supply the values for the first cell.) This solution 
process is referred to as a forward pass, in that we march through the 
computational grid in the direction of increasing (i,j,k). The problem is 
no longer spatially implicit. Each step does, however, require the solu- 
tion of a 5 x 5 system similar to equation (8a) which expresses the inter- 
relation among the Q ' s at each point. Formally, equation (17) represents 
a lower block bidiagonal system which can be solved directly by forward 
substitution, in contrast to the fully implicit block tridiagonal system 
of equation (12). 

The second equation in (14) is treated in exactly the same manner as 
above with the exception that negative eigenvalues are used. This means 
the transport is in the negative (5»n,s) directions and that the dif- 
ferencing proceeds in the negative (i,j,k) directions. For example con- 
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sidering the £ flux through cell (i,j,k), the flux in crosses face (i,j,k) 
but comes from cell (i+1,j,k) while the flux out crosses face (i-1,j,k) 
but comes from cell (i,j,k). This gives 


X i,j,k * 4 ^ A Xl,j,k^l,J,k X i-H,j,k " lT ^ A ^ Q l,j,k^i-1,j.k X i,J,k 
* ^[ c Xj.kk, ,] l.J.k X f.j.k»1 ' 4 ' [C ' (Q L,k»l.J,k-, X l,j.k 


C i, j ,k 


( 18 ) 


or, rearranging 


I 1 - 4 ’! A "Kj,k )] i-,,j.k - ^Xj.k^i.j-i.k- 4 't c '( Q J,j.k)] 1 ( j,k-il 

X l,j,k' X i,J,k' 4T t* ( Q l»l,j,k^l,J,k X i*1,j,k" 4T t B ^ Q i,j»1,k^i,j,k 


x l,J*1.k- 4 ’t C '( tl l.J.k. 1 )]l,j,k X i,J.k*1 


(19) 


In equations (16) through (19) the indices on Q n and X indicate the loca- 
tion at which these quantities are evaluated while the indices appended to 
the brackets denote the faces for which the metrics are to be evaluated. 
Equation (19) represents an upper block bidiagonal system for the unknown 

vector X. . This can be solved directly for the unknowns in a backward 
1 , J » K 
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pass starting at the highest index boundaries and solving in the negative 

p 

(i,j,k) directions since, at each step, the i+1 , j+1 or k+1 values of X 
are known. Again, a 5x5 system similar to equation (8a) must be solved at 
each (i,j,k) location. Appendix B provides a more complete explanation of 
the hierarchy of systems embodied in equations (16) through (19). 

With the solution for X 2 now available, the last of equations (14) 
provides the desired solution vector AQ n , which, in turn, provides the de- 
pendent variable vector at the n+1 time level, 

Q n+1 „ Q n + AQ n (20 ) 
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III. Solution Procedure 

Problem Statement 

The solution process can now be presented. The problem to be solved 
is the sequence of equations in (14) which, when the finite difference 
operations are written out, become equations (16) and (18). As described 
in Appendix B, these equations are actually block bidiagonal systems of 
equations which might be schematically represented (using the first of 
equations (14)) by 



where the subscripts refer to spatial position and n is the total number 
of cells in the computational grid. This system is solved directly and 
quite simply by forward substitution marching through the grid. That is, 

X 1 " X 1 

A 21 X 1 + ( 1+A 22^ X 2 “ “ AtR 2 => ( 1+A 22^ X 2 = * AtR 2 " A 21 X 1 

A 32 X 2 + Cl +A 33 ) x 3 = "AtR 3 => (l + A 33 )X 3 = -AtR 3 - A^ (22) 

A .X + [ 1+A )X = -AtR => (l + A )x = ~AxR - A X 
nn - 1 n-1 v mr n n v nn n n nn 1 n 1 
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But this set of equations is exactly equation (17) (or (19) for the second 
step). Therefore we see that the rearrangement from (16) to (17) or (18) 
to (19) is, in fact, the solution of the block bidiagonal system. The 
value of X in each cell is determined by X in the preceeding cell and R. 

What, then, does the program do? We must realize that each X in (17) 
is actually a vector, X = (x 1 ^.x^ ,xjj ,x^ ) , and each equation in (17) is 
actually a 5x5 system of equations of the form 


1+A 

A 


11 

21 


12 


1+A 
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A._ 

• • • 15 


X 1 


“ R 1i + ^ A 11 x 1 + + A 15 X 5^i-1 



X 2 

= At 

(23) 

• 

• 

1+A 55 


_ X 5_ 


" R 5i + ( A 51 X 1 + +A 55 X 5^i“1 


where the number subscripts now refer to the individual Q's (p ,pu,pv,pw, e) 
and i and i-1 refer to spatial positions. This system must be solved for 
each cell in the computational grid. It is the construction of the indi^ 
vidual terms in this 5x5 system and then the solution of this system, for 
every cell, which is the essential function of the computation. 

Doolittle's Method 

The 5x5 linear systems of equations are solved here using Doolittle's 
method. This method^) is one of a family of techniques in which the 


linear system Ax = b is solved by first factoring the coefficient matrix 
into lower-tri angular and upper-triangular terras, L and U, such that LU - 

■> “I 4 + 

A and LUx = b. This last expression is rearranged to give Ux - L b or Lx 
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= U ^ b which can be solved directly for x by backward or forward substi- 
tution. The aspect of this scheme which is particular to Doolittle's 

method is that the terms on the diagonal of the L matrix are all 1. As 
shown in Appendix B, the 5x5 system involves unknown quantities at the 
cell in question and the known solution from a neighboring cell. The 
backward or forward substitution required to find the unknowns from the 
known solution is a recursive relation which prohibits vectorization (on 
current vector processors) in the direction of the substitution. However, 
a vector solution can be effected by vectorizing along a computational 
line normal to the backward or forward substitution directions. Such 
vectorization has been implemented on a variety of vector processors 
(Cray-IS, Cray X-MP, VPS-32) so that totally vectorized routines for 
solving these linear systems with Doolittle's method are available. 

System Details and Procedure 

In order to properly understand the computational procedure let us 
examine some of the details of the 5x5 system. The first equation in the 
system is 


((1 + A 11 )x 1 + A 12 x 2 + A^x^ + A^Xjj + A 15 X 5 ) i - “AtR^ 

+ (A n x 1 + A 12 x 2 + A 13 x 3 + A^Xjj + A 15 x 5 ) i _ 1 (24) 


Writing out the A's and R's (see equations (17) and (13)) we have 


- *1 £‘ ♦ ♦ !!i) 


3H„ 


Jim 


Jim Jim 


Jim 


(25) 


3Q m 30, 


m 


30 , 


m 


and 

R t n - 5 i F I + 6 i F t + + + 


( 26 ) 
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where i = 1,5,m = 1,5 and the flux vectors (F + , •••, H + ) are functions of 
of (Q., .Q2.Q3.Q14.Q5) * J(p ,pu,pv,pw,e) (see equations (2)). (27) 

Recall that the flux vectors are evaluated at the cell faces; this evalua- 
tion requires that a spatial extrapolation of the Q's from the cell cen- 
ters to the faces be performed. In order to maintain the second order 
spatial accuracy of the overall computation it is necessary to use a 
two-point extrapolation for the Q's which are used to evaluate R on the 
right hand side of the 5x5 system. This extrapolation (see Ref. 1) is 

■ '-SOi.j.k * 

( 28 ) 

■ '-SVlJ.k - °-5«l»2,j,k 

depending on which of equations (14) is being considered. A one point 
extrapolation is sufficient on the left hand side. 

The flux vectors (F,G,H) required in this computation are in the form 
of the subvectors (F ± ,G^,H-). The discussion concerning equation (10) 
describes qualitatively how the flux vectors are split into terms (or 
subvectors) having the eigenvalues as coefficients and how they are 
further divided into the subvectors having positive or negative 
eigenvalues as coefficients. A detailed description of this splitting 
(and the general eigensystem for the Euler equations) is presented in 
Appendix C. From the development of Appendix C the flux vectors are 

K - x’ K,* * x5 K 3 (Cl 6b) 
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t— — 


— — 


— — 

p 


P 


P 

pu 


pu+pci< x 


pu-pck x 

pv 

K - i- 
K 2 2Y 

pv+pcky 

K 3 2Y 

pv-pcky 

pw 


pw+pck z 


pw-pck^ 

pe 

jr-i_ 

■ 

e+p+pc9 k 


e+p-pc0 k 


with 

K = F,G,H when k = C»h,5 

pc 2 + y( (Y - 1 ) e-p<j> ) , <{> = (u 2 +v 2 +w 2 ) , 9 = k u+k v+k w (C7) 

2 K x y z 

A? = 6, , A* 1 = 0 +c|Vk|, Ap = 9 -C|V|, Vk •= (k 2 + k 2 + k 2 ) 1/2 (C9) 

k kk G'l’k c II* xyz 

( ) => divide by|Vk| 


Appendix D describes the construction of the flux Jacobian matrices 
(A,B,C) which are required in this computation. The elements of these 
matrices are formed by a straightforward differentiation of the split 
flux vectors, after the flux vectors have been written as explicit 
functions of the conserved variables, Q. One difficulty associated with 
forming the flux Jacobians is that the flux vectors become long, unwieldy 
expressions when written in terms of the the conserved variables. Another 
tedious aspect of the Jacobian computation is properly assigning 
contributions to the positive and negative components. This assignment 
must be done separately for each cell and each time step. Because of the 
length of the expressions for the flux Jacobian elements, it is not 
appropriate to write them out in this presentation. 
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The time step At is based on the maximum allowable time step in each 
cell volume in order to accelerate convergence for steady state solutions. 
The time step is determined at each point in the grid by (see Ref. 1) 


, At^At^At^ 

At = — - — - 

At^At 11 + At^At^ + At^At 5 

with At* = — Ak - , k = 5,n,E, a - M,5 

max|X k | 


(29) 


To express this in a more convenient form introduce short hand notation 

for the maximum eigenvalues X k = m ^ x 1x^1 and substitute the expression 

m 

for At into At to give 


At = CFL 


AgAnAE 


1 


X X X AgAn + AgAe + AnAg 

^m % ? m X_ X X P X XX, 

E n m K m n m c_ 

mm mm mm 


Simplifying yields 


At = CFL 


A^AnAE 


A£AnX + ACAeX + AnAEX_ 
5 m n m 5 m 


However, in the computational grid A? * An = Ae * 1 so that we have 
CFL 


At 


X. + X + X, 
^m \ ^m 


The eigenvalues, X k> are (see equation (9)) 


\ = k 


H.5 


■k - "x u * V * k z“ “ V V ■ 9 k 1 °l 7k l 

and the maximum absolute value will always be |9 k | + c | V k | so that 


At 


CFL 


I(|M * c|7k|) 


k = 5,n,E 


(30) 
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This is the expression by which the time step is evaluated in the program. 

Returning to the computational procedure we see that the 5x5 system of 
equations at each grid point requires a determination of the flux Jacobian 
matrices (A + ,A - , • • • ,C~) , the time step At and the residuals R n . The 
residuals, in turn, require determination of the two-point extrapolated 
Q's, the flux vectors (F + ,F~, • • • ,H~) , and the sum of the difference opera- 
tions (5^F + +* • •+5| < H - ) . With these quantities determined the system can be 
solved. 

The calculations thus proceed as follows. 

1. two point extrapolated Q's 

2. eigenvalues 

3. flux vectors 

4. sum of flux vector differences 

5. time step 

6. flux Jacobian matrices 

7. coefficient matrix of left hand side 

8. lower/upper decomposition of left hand side coefficient matrix 

9. sum of right hand side terms 

10. solution of the system by forward or backward substitution 

1 1 . update Q' s 
Supporting Calculations 

The above procedure constitutes the essential function of the computer 
program. Additional computations and operations are required, however, to 
obtain a complete solution. These additional steps include establishing 
the initial flow conditions, computing the metrics (that is, the dimen- 
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sions) of each cell, enforcing boundary conditions, determining the vis- 
cous influence and presenting the results in a meaningful manner. Brief 
descriptions of seme of these steps follows. 

The boundary conditions at the far field boundary and the body surface 
are enforced using characteristic variables and phantom points as devel- 
oped in Ref. 1 . In a manner similar to the +/- flux vector splitting this 
method takes advantage of the natural signalling processes (information 
propagation) in the flow to more correctly determine conditions on the 
exterior and interior surfaces of the computational domain. The eigen- 
values are used to determine whether information propagates into or out 
of the computational domain and also to provide expressions for the flow 
variables at the boundaries. Phantom points (fictitous points which are 
immediately exterior to far field boundaries and immediately embedded 
inside solid bodies) facilitate the computations near the boundaries. 
Through extensive numerical experimentation the combination of character- 
istic variable boundary conditions and phantom points appears to provide 
more efficient and accurate evaluation of the boundary conditions than 
other techniques such as extrapolation, zero pressure gradient and use of 
the normal momentum equation. 

The only metrics required in this computational solution to the Euler 
equations are the area vectors of the grid cell faces. These vectors have 
components j(k x ,ky,k z , ) where k ■ S,n,S. For example, Jn z is the z compo- 
nent of area for an n = constant face of a cell. The total area of a cell 
face is J|Vk|; for example, J |Vn| is the total area of an n = constant 
face. The area vectors are computed as the cross product of the diagonals 
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of each face. Also computed are the direction cosines of each component 

of the area vector, which are (k ,k ,k ) = (k ,k ,k )/|Vk|. 

x y z x’ y z I I 

Initial conditions are prescribed as freestream conditions everywhere. 
The flow variables are made non-dimensional using p „„ and c,,,. A perfect 
gas is assumed. The freestream pressure is 

P = p RT 
00 r 00 00 


P » R YR 


P C 
00 00 


which, when made non-dimensional with pc , is 


P« = 


Y 


(31a) 


1 2 

The freestream energy (internal, c^T^, plus kinetic, — q^) per unit volume 

1 , 1 
e oo = y-'j Poo 2 

2 2 2 2 
where q = u + v + w 
^00 00 00 00 

2 2 2 2 
Using to non-dimensionalize and recognizing = q^ /c gives 


e 

00 


1 

(Y-1) 



M 


(31b) 


Numerical Properties 

The scheme as described is implicit, upwind and finite volume. First 
order differencing is used in the left-hand side matrix operator which 


27 



implicitly operates on the dependent variable difference AQ n . This yields 
second-order accuracy in space. For consistency, secondSorder differ- 
encing is used in the residuals on the right-hand side such that the 
overall scheme is second-border accurate in space. 

Because only steady state solutions are of concern here local time 
stepping and first-order temporal differences are used. 

The upwind character of the solution scheme precludes the necessity of 
adding artificial dissipation to damp oscillations as is commonly required 
in central difference schemes. 

Analysis of a scalar equation and a system of equations^ 1 ^ shows that 
the present scheme is conditionally stable. The practical limit for CFL 
number is approximately 20, although under certain flow conditions much 
higher CFL numbers are possible. 
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IV. Viscous Calculations 


Viscous effects are determined through the use of an inverse integral 
method for computation of turbulent compressible boundary layers. The method is 
referred to as integral in that the fundamental basis of the calculation is in- 
tegration of the momentum and kinetic energy equations. Use of a prescribed 
displacement thickness distribution in place of the more usual prescribed pres- 
sure distribution as part of the input to the solution gives rise to the inverse 
nature of this method. The influence of the boundary layers is imposed upon the 
external inviscid flow by a surface source model in which the velocity normal to 
solid surfaces induced by the viscous displacement thickness acts as the effec- 
tive strength of an inviscid source distribution on the solid surface. In the 
present application the viscous calculations are along two-dimensional strips; 
no explicit spanwise calculations are performed. A full description of the pro- 
cedure appears in references (10,11,12). 

Central to this scheme is the representation of the turbulent boundary 
layer by an analytical expression for the velocity profile. This analytical 
expression, developed from curve fits to experimental data, is applicable to 
both attached and mildly separated compressible turbulent boundary layers. From 
this expression the various boundary layer length scales, the skin friction and 
the dissipation can be obtained. The velocity profile used here is actually for 
the equivalent incompressible viscous flow. The compressible flow properties 
are determined from correlations of the three compressible shape factors and 
skin friction with the incompressible shape factor and boundary layer edge Mach 
number. These correlations are also based on curve fits to experimental data. 
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The basic governing equations for the particular inverse integral method 
used here are those for momentum and mean-flow kinetic energy. With suitable 
algebraic manipulations, these two equations become a coupled pair of first or- 
der ordinary differential equations for H (the incompressible shape factor) and 
M 0 (the edge Mach number). The equations contain the four compressible length 
scales, along with M g , c f (the skin friction coefficient), and D (the dissipa- 
tion integral). With the analytic velocity profile, the correlations and the 
displacement thickness distribution (from the previous viscous solution) known 
the pair of differential equations can be solved for the unknowns H and M 0 . In 
order to evaluate the dissipation integral the Cebeci-Smith two-layer eddy vis- 
cosity turbulence model is used for t while the velocity profile provides 
9u/3y. The actual solution of the system is obtained through a fourth order, 
four stage Runge-Kutta routine. 

The solution procedure for a given viscous calculation begins with input of 
H, Ug, Re 0 (the equivalent incompressible shape factor, friction velocity and 
Reynolds number based on 0, the incompressible momentum thickness) to the main 
viscous subroutine (SOURCE). These input quantities come from the previous vis- 
cous cycle. Also input are the conserved dependent variables (p , pu, pv, pw, 
e). A sequence of steps follow which lead to a new or updated velocity profile. 
This profile is then used in the solution of the coupled pair of differential 
equations for H and Mg. Also required in the solution for Hand M g is 5, the 
displacement thickness distribution. The displacement thickness for the current 
viscous calculation is obtained by multiplying the thickness from the previous 
calculation by the ratio of the previous edge velocity to the previous inviscid 
surface velocity magnitude. (On the very first viscous pass, the solution is 
started by specifying the location and displacement thickness at the locations 
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of boundary layer transition and computing a flat plate turbulent boundary 

# — 
layer.) With the 6 distribution available updated values of H and M 0 are de- 
termined. The new values of H and M g distributed along the body can then be 
used with the correlations to obtain distributions of all of the actual compres- 
sible properties of the boundary layer. 

Among the compressible properties determined in the above procedure is the 

£ 

mass flux defect, p e u e 6 , created by the boundary layer. The usefulness of this 
property is apparent if we integrate the differential continuity equation; the 
result is 

(Pv) n « ^ (PeV*) 

where (pv) n is the mass flux per unit area normal to the surface induced by the 
boundary layer and x is the streamwise coordinate. With respect to the external 
inviscid flow, the influence of the boundary layer is to displace the stream- 
lines away from the body. This displacement is equivalent to superimposing a 
mass flux, which is equal to the mass flux defect, normal to the body surface. 

The quantity (pv) n is precisely this normal mass flux. Thus the distribution of 
£ 

PgUgS provides the means to determine the viscous influence on the inviscid 
flow, or, in other words, to determine the viscous-inviscid interaction. The 
normal mass flux applied in the solid wall boundary conditions may be inter- 
preted as a surface source strength. An alternate interpretation of (pv) n is 
that it represents the porosity of the surface. The boundary condition routine 
in this program includes (pv) n in the characteristic variable solid wall compu- 
tations. 
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V. Program Details 


Subroutines 

The solution procedure is directed by a single calling subroutine, STEP. 
Some computations are performed in STEP but its main purpose is to call the 
proper sequence of subroutines to perform the requisite calculations. Perhaps 
the best way to briefly describe each subroutine is to put each into the context 
of its sequence and role in the computational scheme. The following table de- 
tails, in chronological order (in the computer program), each subroutine and its 
function, relevant equation or reference and primary output variables. 


subroutine 

function 

equation 

output variables 

IC 

initial conditions 

31 

p,pu,pv,pw,e 

METRIC 

cell face area vectors 

32 

^X’ 

(k=C,n,5) 

BC 

boundary conditions 

42-47 in Ref. 1 

p ,pu,pv,pw,e 

DELQ 

2 point extrapolated Q's 

28 

p,pu,pv,pw,e 

FLUX 

eigenvalues 

9 

4 (k=5,n,C, 
*=1,4,5) 

FLUX 

flux vectors 

Cl 6, C17 

F^G^H 1 , at 
step n 

DELQ 

residual=sum of flux vector 
differences 

13 

Rj, *=1 .... ,5 

El GEN V 

time step 

30 

At 

FJMAT 

eigenvalues 

9 

\ * 
A k 

FJMAT 

flux Jacobian matrices 

Appendix D 

A ± ,B ± ,C i 

STEP 

coefficient matrix of 5x5 
system. 

17, 19 


AEQLU 

lower/upper decomposition of 
coefficient matrix of left 
hand side 

linear algebra 
text 
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DOOP, DOOM 

right hand side of 5x5 system 

17, 19 

-R n +(A^ .j x. +A 
* * * + ^15 x 5 

DOOP, DOOM 

forward/backward substitution 
solution of 5x5 system 

linear algebra 
text 

AQ n 

STEP 

update Q's 

20 

p ,pu,p v,pw, e 

BC 

boundary conditions 

H2-47 in Ref. 1 

p ,pu,p v,pw, e 

SOURCE 

viscous calculation 

Section IV 

pv n 

PVAR 

print results 


Cp , x/ c , Cj_ , c 

Variables 





As would be expected in a program designed to solve a complex system of 
equations, there are a large number of variables. The names used for the vari- 
ables, in most cases, are representative of the variables' physical or mathe- 
matical meanings. Brief definitions of the more important quantities are pre- 
sented here. 


The following variables appear in nearly all subroutines. 

I,J,K = C,n,?; indices for grid points 
NI ,NJ ,NK « number of (I,J,K) lines 

R,RU,RV,RW,E,P = p ,pu,p v,pw, e,p in Cartesian frame, (x,y,z) 

X,Y,Z = grid point coordinates in Cartesian frame, (x,y,z) 

AIX,AIY,AIZ - Cartesian components of the area vector of an I constant cell face 

AJX,AJY,AJZ « Cartesian components of the area vector of a J constant cell face 

AKX,AKY,AKZ = Cartesian components of the area vector of a K constant cell face 

SADAI ,SADAJ ,SADAK = magnitude of the area vectors for I,J, or K constant cell 

faces 


The following variables appear in only some of the subroutines: 

B (L ,M,N ,1 , J ,K) = flux Jacobian matrices; L=1 ,6 and refers to A + ,...,C ; M=1,5 
and refers to the flux vector component; N=1 ,5 and refers to 
the dependent conserved variable; I,J,K refer to grid location 
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D(L,M,N,I,J,K) = coefficient matrix on the left hand side of the 5x5 systems; 

L=1 ,2 and refers to the forward or backward passes; M=1,5 and 
N=1 ,5 refer to the location within the matrix; I,J,K refer 
to grid location 

DR , DRU , DRV , DRW , DE = right hand side of two pass algorithm (R.X 1 ); AQ n after two 

pass algorithm 

DT = time step 
In DELQ, 

RR,RUR, . . . ,RWL,EL » right and left, two point, extrapolated dependent variables 
XR ,XRU,XRV ,XRW,XE = flux vectors 

In DOOP and DOOM, 

T = sum of transport terms on right hand side of the 5x5 systems 

Z - total right hand side of the 5x5 systems, initially; solution to 5x5 system, 
finally 

Y - intermediate variables in solution of 5x5 system 
In FJMAT , 

EV1,EV4,EV5 * eigenvalues 

CG1,CG2,CG3 «• switches to assign terms to appropriate (±) flux Jacobians 
AXT,AYT,AZT - k x ,k y ,k z 
In FLUX, 

EV1R, . . . ,EV5L - eigenvalues predicted on cell faces using right and left extra- 
polated dependent variables 

XR , XRU , XRV , XRW , XE - flux vectors 

In EIGENV , 

CONU , CONV , CONW = 9 k 
El ,EJ ,EK = | Q k | + c | Vk | 

DT = time step 
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Input parameters (in order of read), 

CFL ■ CFL number 

FSMACH = freestream Mach number 

ALPHA, BETA, PHI = roll, pitch (angle of attack) and yaw angles 
NB * number of printouts 

NT = number of computational cycles per printout 

NV * number of computational cycles per call to viscous calculations 
ITL,ITU,ILE = I values at lower and upper trailing edges and at leading edge 
KTIP - K value at wing tip 

XI = array storing the computational grid coordinates 

IFREQ » frequency of calls to EIGENV, FJMAT and AEQLU (not read but set in MAIN) 
NI,NJ,NK = grid size (not read but set in PARAMETER statements) 


Output, every 5 cycles 

NCYC = cycle number 

RTMAX - maximum DR 

RTRMS = rms value of DR 

ETMAX * maximum DE 

ETRMS = rms value of DE 

XL2 - sum of squares of DR,...,DE 

TCL = lift coefficient (wing) 

TCD = drag coefficient (wing) 

NSUP = number of supersonic points 

Output, every NT cycles 

ZLOC = spanwise position 

S = chordwise position 

CP2 = surface pressure coefficient 

CL *« sectional lift coefficient 

CD = sectional drag coefficient 

TCL = wing lift coefficient 

TCD = wing drag coefficient 
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Execution sequence 


MAIN 

call IC 
read input 
call METRIC 
call BC 
call PGEOM 

do 3 L=1 ,NB 

do 2 M=1 ,NT 

call STEP 
call DELQ 

call FLUX, (NK-1)*(NJ-1)+(NK-1)*(NI-1) + (NJ-1 )* (NI-1 
call EIGENV, NB*NT/IFREQ times 
call FJMAT, 6*NB*NT/IFREQ times 
call AEQLU , 2*NB*NT/IFREQ times 
call DOOP 
call DOOM 
call BC 

call SOURCE NB*NT/NV times 
call INVBL 
call CORREL 
call LININT 
call RUNGE 

2 continue 
call PVAR 

3 continue 


times 
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VI. Applications 


This code has been used to compute transonic flow past a variety of 
two-dimensional and three-dimensional bodies. To demonstrate the applica- 
tion of the program, results from computations of a representative wide 
body subsonic jet transport will be presented. 

Geometry and Grid 

The configuration in question is a complex wing/fuselage combination 
designated for this work as the Pathfinder geometry. The wing is a high 
aspect ratio, compound sweep surface with supercritical sections, span 
wise washout and positive dihedral. This wing is mid-mounted on a simple 
axisymmetric fuselage which has a slightly blunted ogive nose and a blunt 
base without boattailing. For the computational geometry there is no 
wing root fillet or fairing. A geometry very similar to this has under- 
gone wind tunnel testing at NASA Langley and some results from those 
experiments will be used for comparison. 

To generate the computational grid the program FL059, authored by T. 
Jameson, ^3) was applied. This code is a full Euler equation solver but 
also includes a grid generation package for wing/body/ tail geometries, 
which has been slightly modified by Dr. J. Luckring of NASA Langley and 
further modified during the present work. Figures 1 through 6 show por- 
tions of the grid produced by FL059 for the Pathfinder configuration. 
This is a C-H grid with 97 I lines and 17 each of the J and K lines. The 
lower and upper wing trailing edges are at I = 19 and 79 respectively; the 
leading edge is at I = 49 . The wing tip is at K - 11. 
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The K = 1 surface is the plane of symmetry (z = 0) above and below 
the body and lies on the surface of the fuselage, cutting through the wing 
at the wing/fuselage junction. Figure 1 shows the K = 1 surface (viewed 
along the z axis) in the immediate vicinity of the fuselage. The blunt 
base, wing root profile and fuselage nose shape are clearly evident. An 
oblique view of the K * 1 surface near the fuselage appears in Figure 2; 
prominent in this figure is the rapid transition at the fuselage base of 
lines on the surface to lines on the mid plane. On Figure 3 an enlarge- 
ment of the wing root area (K = 1 , view along z axis) is presented. The 
supercritical section is clear. Two less desirable grid properties are 
also apparent. There is a slight irregularity in the grid lines immedi- 
ately downstream, of the trailing edge. More importantly, the cells on 
the surface (upper and lower) on the aft half of the airfoil are somewhat 
large; this may adversely influence computations on the compression (down- 
stream facing) surfaces of the wing. 

Three views of the J = 1 surface are presented in Figures 4, 5 and 6. 
The entire J • 1 surface appears in an oblique view in Figure 4 where the 
left-most line runs along the wing leading edge. The grid extends ap- 
proximately one and one-half semispans beyond the wing tip in the z direc- 
tion and approximately one semispan downstream of the fuselage base. The 
abrupt bend in the bottom most grid line is at the fuselage base. In 
Figure 5 the wing is viewed along the y axis. Here the slight change in 
leading edge sweep is apparent. Another oblique view of the wing, but 
somewhat enlarged, is presented in Figure 6. The geometry of the trailing 
edge can be more clearly discerned in this view. 
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The outer boundaries of the grid are relatively near the body, being 
roughly one semispan in each direction from the body, (The extent in the 
y direction and ahead of the body are not shown but are about one and 
one-half semispans each.) Numerical experiments have not shown any obvi- 
ous adverse influence resulting from insufficient distance to the outer 
boundaries although this spacing and its effects may need to be further 
explored. 

Calculations 

Computations were performed over a range of values for the various 
parameters which can be set in the code, in particular IFREQ, CFL, NV , 
NB*NT, DST, XTOP and XBOT. The best combination to minimize machine time, 
maximize convergence rate and still obtain reasonable output involved 
running for 100 to 300 cycles (NB*NT) with IFREQ=9999, CFL=15 and NV 
between 10 and 20. This means only one call is made to FJMAT, AEQLU and 
El GEN V and from five to 30 calls to the viscous routines (when they are 
used) . For this combination the residuals are reduced by at least three 
orders of magnitude and the lift is within five to ten percent of its 
asymptotic value. All results presented here are for parameters in these 

ranges. The examples to be presented here are for two cases with the 
Pathfinder geometry: case (1), M = 0.82, a = 2°, Re_ = 9.9*1 0^ and case 

oo C 

(2), M = 0.70, a = 2°, Re_ = 5.3*1 0^ . These cases correspond to wind 

oo C 

tunnel experiments conducted at NASA Langley Research Center. 
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Figure 7a presents computed pressure distributions for case (1 ) at 
the 45? span location along with an experimental distribution at the 43.2? 
span location position. The two curves on this figure are for an inviscid 
run of 100 cycles and a viscous run of 100 cycles with five calls to the 
viscous routines. For the viscous calculations the turbulent boundary 
layer was begun at x/c values of 0.05 and 0.2 on the lower and upper 
surfaces, respectively, while the initial displacement thickness (6*/c) was 
0.0010 and 0.0012, lower and upper. On Figure 7a, the lower surface 
agreement is quite good up to x/c = 0.8; the deviation downstream of that 
point is expected in view of the relatively coarse grid. On the upper 
surface, the leading edge suction peak is properly computed and the basic 
character of the recompression is represented by the calculations. The 
shock, however, is not accurately captured. This may be a consequence of 
the coarse grid. 

Distributions at other spanwise locations for case (1) appear in 
Figure 7b-f. The agreement between experiment and calculation deterio- 
rates as the root and tip are approached. The nature of the distributions 
outboard of the 45? location suggests that the wing twist is not correctly 
modelled. There is progressively too much upper surface expansion and too 
little lower surface expansion as the tip is approached which indicate 
that these sections are at too high an angle of attack. At the root there 
is insufficient upper surface expansion indicating that the angle of 
attack is not large enough. It is most likely that angle of attack is not 
the only problem at the root. The interaction between the fuselage and 
wing is probably not being adequately computed also. 
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As an example of the influence of angle of attack, Figure 8 is pre- 
sented in which viscous runs (same parameter values as in Figure 7) for 
three angles of attack are compared. It should first be observed that the 
experimental angle of attack is 2° but that the results in Figure 7 are 
actually at 2.5° as this gives reasonably good agreement at the 45? span 
position. Figure 8 shows computations at 2°, 2.3° and 2.6°. At the 
outboard sections, the tip in particular, the calculations more closely 
match the experimental results as the angle of attack is decreased. The 
trend is reversed at the inboard sections. These results further support 
the idea that the wing twist is not correct. 

Results for case (2) appear in Figure 9, where the computations are 
viscous and were run for 250 cycles with NV=10. The agreement here is 
good except at the most inboard section. A spanwise variation which is 
probably due to incorrect twist is still detectable but, overall, the 
calculations quite satisfactorily predict the flow. 

Concluding Remarks 

Development and application of a code to accurately compute 
three-dimensional, transonic flow with viscous effects is not a simple 
matter as the preceding discussions hopefully show. The emphasis during 
this research has been to concentrate on refinement of the numerical 
scheme with particular attention paid to improving the computational speed 
and eliminating numerical oscillations. Both of these goals have been 
closely approached, if not fully met. However, there is clearly much that 
still needs to be addressed. An improved and finer mesh is obviously 
necessary to obtain accurate flow predictions. The grid size chosen here 
was dictated primarily by economy rather than accuracy. A fully three- 
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dimensional viscous calculation scheme should be implemented in order to 
accurately compute root and tip flows. Multi-grid and block schemes and 
adaptive grids should also be considered. 

Although many improvements can be made, the code as described here is 
an efficient and robust Euler equation solver, relatively easy to use and 
with reasonable accuracy with a coarse mesh. When coupled with a well- 
developed computational grid this program, in its present form, should 
prove to be a useful tool for engineering calculations of three-dimen- 
sional, high Reynolds number, transonic flow. 
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Appendix A. Coordinate Transformation 


The steps required to transform the Euler equations from Cartesian to 
curvilinear coordinates are described here^. The logic behind this proce- 
dure may be more clear if it is realized that the main purpose of this 
transformation is to replace all Cartesian derivatives with curvilinear 
derivatives. 

The Euler equations in Cartesian coordinates are in the form 

q t + f x + g y + h z " 0 (A1) 

Consider ( 5 , n, c) = function (x, y, z, t), 

(x, y, z) = function ( 5 , n, S, x) and t = t. 

Cartesian derivatives are then 


9_ 

at 

a 

9x 


3 a 9 3_ 

T t 9x 5 t 95 n t 9n + 5 t 9c 


T x 9x + ^x 9£ + n x 9n ^x 9x ’ etG< 


(A3) 


and curvilinear derivatives are 

9 v 9 ^ 9 , 9 9 

9x = ^x 9t X C 9x y x 9y Z x 9z 

9 .9 9 9 9 

- t_ -T- + X_ -r- + y_ -r— + Z_ -5- , etc.* 

9^ 5 9x 5 9x ■'S 9y 5 9z 

Substitute (A3) into (A1), realizing t * x * f(x, y, z) gives 


(Ait) 


q T * 5 t Q 5 * V„ * 5 t q t * * C x f t 

* E y 8 S * n y s n * 5 y h ; * { A * "zV Vs ‘ 0 


(*5) 


H5 



At this point we have curvilinear derivatives of the dependent vari- 
ables and flux terms (q, f, g, h) but we still have Cartesian derivatives of 
the independent variables (£, n, c, t) . We must, therefore, find expres- 
sions for the Cartesian derivatives in terms of curvilinear derivatives. To 
do this first write equation (4) in matrix form 

[curv deriv] ■ [j] [cart deriv] (A6) 

and realize that the determinant of the coefficient matrix, | [ J ] | is the 
matrix Jacobian, here denoted simply as J, 


x ? ( Vc ' ‘ y 5 ( Vc ' Vt> * 




(A7) 


Multiplying (A6) by [J] 1 gives the desired relation 


[cart deriv] » [j] -1 [curv deriv] 


(A8) 


Apply (A8) to £ t , £ x , etc. to give 
£t = -x t^x "* ^T^y ~ z t^z 

(A9) 

c x - J_1 (y n z ? ■ Vc 5 * and so on * 

Before equations (A9) are used in equation (A5), a simplification is 
possible. This simplification is achieved by multiplying (A5) by J and 
using the product rule for derivatives in the following particular form. 
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(A10 ) 


Jk n, i ■ h (J V) - c Ik 

where k = (5, n, ?, t), m = (x, y, z, t) and C = (q, f, g, h) , respectively. 
Doing so yields an equation of the form 

If [ J o]' * 1 1 E J C«t<5 * V * 5 y g * «z h)1 2 * 3 * 4 

- 1 & (•» * It 4 If < J V * ij 5 (1,,) 

- f Ife (J V 4 1 ; <J \> 4 lr W! x>' 6 * 7 4 8 

Terms 3 and 4 are just like 2 but with n or c replacing 5. Similarly, terms 
7 and 8 are like 6 both with g and y or h and z replacing f and x. 

The bracketed quantities in terms 5 through 8 of this equation are the, 
so-called, invariants of the transformation and are zero. This can be shown in 
a 'brute force' sort of way as follows. Consider the bracketed quantity in 
term 6: 


k (JE x> 4 k (J V 4 k (J! x> 


chain rule 


_ 3 „ 3J . 3 3J T 3 3J 

J 35 C x + 5 x 35 + J 3n n x + n x 3n J 3? ? x ? x 3? 


rearrange 


r 3 3 3 i _ 3 J 3 J 3 J 

= J ^35 5 x + 3n n x 3c 5 x^ 5 x 35 n x 3n C x 3c 

change order of differentiation and use (A3) 


J l;(i^ 4 §;" 4 i^> 4 -H J 
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evaluate 


■ J h <3) * h ( Wt ■ Vn’ * "• ) 

- 0 + 0 (if change order of differentiation). 

Terms 5, 7 and 8 behave similarly. 

Equation (All) is now of the form 

|Q + |L + |2 + |h .o 

3t 3^ 3n 3? 

where Q = Jq, F « j(5 t Q ♦ 5 x f + 5 y g + £ z h) 

G = J(n t q + n x f + n y g + n z h), H ■ J(c t q + ; x f + c y g + C z h) 


(A 1 2 ) 


(A13) 


The quantities 5 f> £ x , •••, c y » t z can be replaced by equations (A9) to 
finally yield an expression entirely in terms of curvilinear derivatives, 
and the transformation is thus complete. 


Appendix B. Equation Hierarchy 


Equations (16) and (18) actually represent systems of equations imbedded 
in a system of equations. That is what is meant, for example, by the descrip- 
tion that (16) is a lower block bidiagonal system. The block is, itself a 
system of equations. To show this clearly, let us consider a different fac- 
torization of equation (12) which will yield the same hierarchy of systems as 
in (16) and (18) but with fewer terms. The structure of the systems will then 
not be obscured with lengthy expressions. 

Equation (12) can be factored as 

(i + Ax6 .A* )(l + Ax5.A”)(l + At< 5 .B • ) (i + A5 B~)(l + 6 C-)(l + 6 C~)AQ n = 

1 1 j J K K 

-AxR n (B1 ) 

and solved in a six step scheme 

(I + AxA^tjx 1 " -AxR n 

(i + Ax6.A~)x 2 - X 1 


(i + Ax6 k C~)x 6 = X 5 


AQ n = X 6 


(B2) 


The motivation for this factorization and scheme and the solution process are 
exactly the same as described for equations (13) and (14). 

Let us consider the first equation in (B21), 

(i + Ax6 i A»)x 1 = X 1 + Ax6 j ,(A + X 1 ) = -AxR U 03) 
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and write out the difference operator with indices as described for equations 
(16) - (19). 


X + 

i.j.k 


At ^ A ( Q i,j,k^i,j,k X i,j,k " At ^ A ^ Q i-1 , j ,k^i-1 , j ,k X i-1 , j ,k 


-AxH 


n 

i,j,k 


(B4) 


We must realize several facts concerning equation (B4) and its notation. The 
(i,j,k) subscripts refer to spatial location, the 1 superscript to the solu- 
tion pass and the n superscript to the time level. In addition, at each 
(i,j,k) location there are 5 distinct Q n (representing p ,pu,p v,pw, e) and 

1 + + 
consequently 5 distinct X . Finally, A denotes so that there are 25 A 

+ 

terms since there are 5 flux vectors F^. 

To display the structure clearly, suppose there are only 6 cells in the i 
direction. Also suppose we are evaluating (B4) along a line of constant j and 
k so that the (j,k) indices can be dropped. Further, let us drop the super- 
scripts for this example and finally let A stand for AxA + . Then (B4) becomes 

X. + A.X. - A. ,X. = - AxR. with i = 2,6 (B5) 

l ii 1-1 i-1 i 

Writing (B5) out gives 

X 1 " X 1 
A 1 X 1 + (l + A 2 )X 2 = -AxR 2 

- a 2 x 2 + (i + a 3 )x 3 = -AxR 3 

- a 3 x 3 + (i + a 4 )x 4 = -AxR 4 

- a 4 x 4 + (1 + a 5 )x 5 = -AxR 5 
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- a 5 X 5 * (' * V X 6 ■ ‘ a,R 6 

or, in matrix form 









“ “ 


i ' 

> X 

1 

0 

0 

0 

0 

0 


X 1 


A 1 

1+A 2 

0 

0 

0 

0 


X 2 


R 2 

0 

A 2 


0 

0 

0 


X 3 

= -Ax 

R 3 

0 

0 

# 3 

1 + Aj| 

0 

0 


X 4 


R 4 

0 

0 

0 

A 4 

1+A 5 

0 


X 5 


R 5 

0 

0 

0 

0 

A 5 

1+A 6_ 


X 6_ 


R 6 


(B6b) 


The lower bidiagonal form of (B6b) is now clear. Remember that the subscripts 
1 through 6 here refer to spatial locations. 

Now let us consider the details of any one of the equations in (B6a) or 
(B6b). That is, we now have a fixed value of i and therefore are considering 
a particular location in space. Since (B5) represents a typical equation in 
(B6), we can use (B5) realizing that i is fixed and rearrange to give 


(l + A i )X i 


-AtRi ♦ A 1 - 1 X 1 . 1 


(B7) 


But, at each i there are 5 distinct components of X (say, x^, t = 1,5) and 25 
distinct A terms. Thus (B7) actually expresses the following matrix equation 


1+A 11 A 12 A 1 3 A 14 


A 21 1+A 22 



15 


1+A 


55 4 


| R 1 + 1 X 1 +A 12 X 2 +A 13 X 3 + A l4 X i| +A 15 X 5^i~1 




(B8) 


-R 5 . + ^ A 51 X 1 + A 52 X 2 +A 53 X 3 +A 5 1 < X 4 +A 55 X 5 ^i-1 
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Here the number subscripts refer to the dependent variables ( = 
(p ,p u,p v, pw, e) while the i subscripts refer to the location, which is now 
fixed. So, each equation in the factored scheme (B2) represents a system of 
equations (B6b) over space. Each equation in system (B6b), in turn, repre- 

sents a system of equations (B8) over the dependent variables. This last 
system (38) is the "block" in the block bidiagonal descriptor. Notice that 

this system (B8) is similar to system (8) in the main text. 

The procedure described above would be carried out for each of the six 
steps in system (B2), with appropriate care in evaluating the difference 
operators as described for equations (17) and (19). As a note, the six step 
scheme of (B2) satisfactorily provides solutions to the Euler equations. It 
requires more computation than the two step scheme in the present program but 
requires less memory and thus is advantageous when memory limitations are 
significant. 

Relating these ideas back to the present scheme, we see that equations 
(14) correspond to (B2), equation (16) to (B4), and (17) to (B7 ) . Thus, the 
first equation in the factored scheme (14) represents a lower block bidiagonal 
system of equations, (16). Each equation in system (16), in turn, represents 
a system of equations (17) over the dependent variables for fixed (i,j,k). In 

a similar fashion, the second equation in (14) represents an upper block 

bidiagonal system which is expressed in equation (18). For fixed (i,j,k) each 
equation in system (18) represents a system of equations, represented by 
equation (19), over the dependent variables. 
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Appendix C. Eigensystems of the Transformed Euler Equations 


The development of the eigenvalues, eigenvectors and split flux vec- 
tors of the transformed Euler equations is briefly presented here. This presen- 
tation will be in condensed notation to emphasize the essential steps. Details 
of the development, particularly in regard to the actual elements in the various 
matrices are available in References 1 and 2. In addition, several theorems and 
operations of matrices and linear algebra are required and these are documented 
in any quality text on linear algebra or numerical analysis. 

To aid in the description it is helpful to first define some notation. 
Certain terms represent any one of three terms depending on the direction (in 
computational space) of interest. These are 

k = 5,n,£ = curvilinear directions 

K = F , G,H = flux vectors (see equation (2) in main text) 

K = — = A,B,C = flux Jacobian matrices, where the Q are the dependent 
variables in the conservative form of the Euler equa- 
tions (eqn (2)) 

k = a,b,c = coefficient matrices for the nonconservative Euler equations 

Typically these are used such that when k « 5 then K=F , K=A and <=a, and so on 
for k=n and k=c. Another item of notation here is that when (t,5,ri,c) appear as 
subscripts they refer to derivatives (3/3t, 3/35. 3/3n, 3/3c). similarly, 

(x, y, z) subscripts mean (3/3x, 3/3y, 3/3z). All other subscripts will not indi- 
cate differentiation. 

We may now proceed with the development. The starting place is the Euler 
equations in conservative form and in curvilinear coordinates, 

Q t + F C + G n + H C = 0 (C1 } 

This can be written in the quasilinear form 
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0 


(C2) 


Q t + AQ^ + BQ n + CQ 5 - 
where, recall, A = 9F/3Q, B = 3G/3Q, C ■ 3H/3Q are terms of K. 

With the governing equation written as (C2) the eigenvalues of the govern- 
ing equation are, in fact, the eigenvalues of the flux Jacobian matrices K. How- 
ever, when K is actually written out (see Appendix D) it is seen that there are 
very few zero elements in K and consequently that it would be a difficult task 
to extract the eigenvalues of K. To find the eigenvalues, a similarity trans- 
formation (that is, one which preserves the eigenvalues) of K is sought which 
yields a matrix with many more zero elements and hence makes extraction of the 
eigenvalues of Kmore tractable. This transformation is achieved by considering 
the nonconservative form of the Euler equations 

q T + aq^ + bq n + cq^ = 0 (C3) 

where q - j[p ,u, v, w, p] T . 

Now (C2) can be written as 

Mq T + AMq^ + BMq^ + CMq ? - 0 (C4) 

where M « 3Q/3q. 

Multiply (C4) by M -1 to give (with I ■ identity matrix) 

Iq T + M _1 AMq^ + M i,1 BMq n + M -1 CMq c = 0 (C5) 

Comparing (C3) and (C5) we see that 

a = M _1 AM, b *= M~^BM, c - M _1 CM (C6a) 

and, therefore, ic * M^KM, K *» MkM - 1 (C6b) 

Equation (C6b) is the required similarity transformation. The matrices K and < 
are similar and therefore have the same eigenvalues. Writing out k we have 
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where 9 k =k x u+k y v+k z w 
pc 2 = T C (Y-1 ) e-p<t« ) 

(u 2 +v 2 +w 2 ) (C7) 

The eigenvalues, X k , are found from the solution of 

|k - IX k | - 0 (G8) 

and are 

A k " *k " X k " 6 k’ *£ “ 0 k + c l Vk !« X l " e k - c l Vk l (° 9) 

with | Vk | - (k 2 + k 2 + k|) 1/2 . 

The eigenvalues have been found and now the task is to find the eigen- 
vectors of K. The eigenvectors, T k , corresponding to the eigenvalues are such 
that 

K - VJ" 1 (CIO) 

where A k is the diagonal matrix whose diagonal elements are the eigenvalues, 
(C9). To determine T k , consider a similar expression for tc 

< - PkVk 1 (C11) 

The matrix P k is also composed of eigenvectors corresponding to the eigenvalues 
(C9). Here the columns of P k are the right eigenvectors of k, each column being 
a linearly independent set of eigenvectors. These eigenvectors are normalized 
such that multiplying the matrix of the right eigenvectors times the matrix of 
the left eigenvectors gives the identity matrix. (As explanation, given a ma- 
trix [A] and [A]x r = Xx p then X is an eigenvalue and x p is the corresponding 
right eigenvector. Similarly given [A] and Xj_[A]=Xj > X then X, again, is an eigen- 


0 

0 

0 

0 


P k > 

e k 

0 

0 


K x pc £ 


P k y P k z 

0 0 

6 k 0 

o e k 

kyp c 2 k z p c 2 


x/p 


Vp 


c z/p 
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value and x^ is the corresponding left eigenvector.) It is much easier to find 
P k (refer to a good linear algebra text), than T k . Once P k is found then T k can 
be readily obtained by: 

from (C 6 b) K = MkM -1 and using (Cl 1 ) then, K = MP (< A k P k 1 M 1 

so that T k = MP k , T " 1 = M^P " 1 (Cl 2) 

The eigenvalues and eigenvectors of the flux Jacobian matrices are now 
available in (C9) and (Cl 2). It remains to obtain the flux vectors as a linear 
expansion in the eigenvalues (eqn (10)). Because the flux vectors are homogene- 
ous of degree one in the Q's we can write, by invoking Euler's Theorem 

K = KQ (Cl 3) 

(Essentially, a homogeneous function f(x), say, of degree n is one such that 
f( ax) = a n f(x) where a is a constant. Euler's Theorem essentially says that for 
f(x) continuously differentiable and homogeneous of degree n, then x df/dx « nf. 
These ideas can be extended to functions of several variables.) 

Using (CIO) in (Cl 3 ) gives 

K = T k A k T?Q (Cl 4) 

We can write the diagonal matrix A k as the sum 

A k " A k I 1,2,3 + X k X 4 + x k I 5 (C15) 

where I ^ 23 * 3 a matrix with 1 as the first 3 diagonal elements and 0 

elsewhere, 1 ^ has 1 as the fourth diagonal element with 0 elsewhere and 

I 5 has 1 as the fifth diagonal element with 0 elsewhere. Then 

K ■ 4Vl,2,3 T ‘k' a * ‘kW?® * »kWi;' 0 (c,6a) 

= A k K 1 + xJk 2 + X k K 3 (Cl 6 b) 
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Equation (Cl 6) is the expansion of the flux vectors in the eigenvectors and is 
the desired split form. When the details of K 1 .Kg and are determined we fi- 
nally have 


K 1 = J 


Y-1 


~p ~ 


P 


P 

pu 

J 

pu+pck x 

J 

pu-pck x 

pv 

"2 = 27 

p v+ p c k y 

K 3 = 2Y 

pv-pck y 

pw 

p<f> 


pw+pck z 


pw-pck z 

e+p-pc9 k 

(Y*1)_ 




(Cl 7 ) 


where (k x ,k ,k z ,9 k ) = (k x ,k ,k z ,9 k )/|Vk| and <J> and 9 k are defined with equation 


(C7) 



Appendix D. Evaluation of the Flux Jacobian Matrices 


The dependent variables of the conservative, transformed Euler equations 

are 

Q = Q m , m = 1,5 * j[p,pu,pv,pw,e] T = [ Jp ,Jpu,Jpv, Jpw, Je] T (2) 

The flux vectors are K = F,G,H and are given by 

K = A J C K 1 + A JJk 2 + XjjK 3 (Cl 6 b) 

where the subvectors K 1 2 3 are given in equation (C17) and are given by 

equation (C8). The flux Jacobian matrices are formed by differentiating K with 
respect to Q, 

A lm ■ 3 V 3 «m- B,„ - C - a^/aa,, (6b) 

In order to carry out the differentiations, it is necessary to write the flux 
vectors explicitly in terms of the conserved variables, Q. For example, the 
second element in the K 2 subvector, K ?5 , is 

K 22 3 fy (P u + PCk x ) = iy( J P u + J P°k x ) (C17) 

Now, from (C7) 

pc 2 = Y((Y-1)e - p<j>) = Y(Y-1)(e - J p(u 2 + v 2 + w 2 )) 
and we can write 

Jpc = (Jp*Jpc 2 ) 1/2 = (Y(Y-1 )) 1/2 (Jp * J(e - p(u 2 +v 2 +w 2 )) 1/2 (Dla) 


or 
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Jpc - (Y(Y-1)) 1 / 2 (Jp*Je •- -1 ((Jpu ) 2 + ( Jp v ) 2 + (Jpw ) 2 )) 1/2 (DID) 

*= (Y(Y- 1 )) 1 / 2 (Q 1 Q 5 -|(Q§ + Q 3 + qJ )) 1/2 (using (2)) (Die) 

So, K 22 becomes 

K 22 “ ^y (Q 2 + (T(Y-1)) 1 / 2 (Q 1 Q 5 - ~(Q| + Q 3 + Q ^)) 1/2 k x ) (D 2 ) 

4 

We now turn our attention to X k , 

X k - 0 k + c|Vk| • k x u + kyV + k 2 w + |Vk| p ~ (D3a) 

= (k x Jpu + kyJpv + k z Jpw + |Vk|Jpc) (D3b) 

Using ( 2 ) and (Die) in (D3b) gives 

X k - ^-(k x Q 2 +k y Q 3 +k z Q 4 +|Vk|(Y(Y-1)) 1 / 2 (Q 1 Q 5 - ^ (qI+Qj+qJ) ) 1 /2 ) (D4) 


Finally, to construct the contribution of K 22 to K we must multiply (D2) and 
(D4) together. Clearly, the contribution is a lengthy expression. Moreover, 
this product is only one of the fifteen terms which compose K from (Cl 6 b). Suf- 
fice to say, the total expression for K is quite involved. 

Once all fifteen terms of K are constructed, in the manner described above, 
the derivatives with respect to the five Q’s can be obtained. Now, we must re- 
member that the solution algorithm is based not only on splitting the flux vec- 
tors into terms with eigenvalues as coefficients but the fluxes are further 
split according to the 3 ign of the eigenvalue, as expressed in (10c). The flux 
Jacobians are correspondingly split, as in (11). Consequently the fifteen terms 
of K must be separated into terms with positive and negative eigenvectors; dif- 
ferentiation then yields the appropriate flux Jacobian matrices A*, B~ and C*. 
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In the program the derivatives of all fifteen terms for each of A, B and C 
are computed. Coefficients of either 1 or 0 are then assigned to each deriva- 
tive according to whether the corresponding eigenvalue is positive or negative 
and which flux Jacobian, the plus or minus, is being evaluated. This procedure 
is organized in an efficient and orderly, although quite lengthy, scheme. 
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Figure 1. Portion of the K-l siirface for the Pathfinder geometry. 





Figure 3. Closeup of the root section of the K=1 surface. 






Figure 7a. Section pressure distribution, ^=*0.82, a=2°, Re=9.9*10 
E - experiment, 43.2% span; I = inviscid, 45% span; 

V «* viscous, 45% span. 
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Figure 7f. E, 96.1%; I and V, 95% 



Figure 3a. Section pressure distribution, Hx>=0.82, a=2°, Re=9.9*10 6 ; 

E = experiment, 13.1% span; (l)=viscous, 2°, (2)=viscous, 
2.3°, (3)=viscous, 2.6°; all at 15% span. 
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Figure 8b 
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Figure 8c. 
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Figure 8e. 









Figure 9a. Section pressure distribution, M»=0.7, a=2°, Re=5.3*10 ; 

E = experiment, 13.1% span; , viscous, 15% span 
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